function fval = int_prod_disp(L, B, a_b, a_s, EPS_L_BOUND, EPS_H_BOUND, S_PARAM, M_PARAM, THETA, DELTA, CHI, tau_k, weight)

sig = S_PARAM;
mu  = M_PARAM;
f   = @(x, a) 1 ./ (x * sig * sqrt(2 * 3.14159265359))...
    .* exp(- (log(x) - mu - a).^2 / 2 / sig^2)...
    .* (x >= EPS_L_BOUND) .* (x <= EPS_H_BOUND);
int_method = 'auto';
Rel_Tol    = 1e-5;
Abs_Tol    = 1e-6;

F   = @(x,a) logncdf(x, mu - a, sig) .* (x >=EPS_L_BOUND);

% lower bound
fun_l_bound      = @(a) EPS_L_BOUND;

% upper bound
fun_h_bound      = @(a) EPS_H_BOUND;


fun_M = @(x, L, a)F(x, a) - integral(@(eps_t)...
    (F(max(L / (1 - THETA) -  THETA / (1 - THETA) * eps_t, eps_t), a) - F(eps_t, a)) .* f(eps_t, a),...
    fun_l_bound(a), x);
fun_m = @(x, L, a) f(x, a) - f(x, a) .* (F(max(L / (1 - THETA) -  THETA / (1 - THETA) * x, x), a) - F(x, a));
fun_g = @(x, L, a) fun_m(x, L, a) / fun_M(fun_h_bound(a), L, a);

moment_1   = (1 - weight) * exp(mu + sig^2 / 2) +...
    weight * integral(@(x) x .* fun_g(x, L, a_b), fun_l_bound(a_s), fun_h_bound(a_b));

moment_2   = (1 - weight) * ((exp(sig^2) - 1) * exp(2 * mu + sig^2) + (exp(mu + sig^2 / 2))^2)+ ...
    weight * integral(@(x) x.^2 .* fun_g(x, L, a_b), fun_l_bound(a_b), fun_h_bound(a_s));

fval   = sqrt(moment_2 - moment_1^2) / moment_1;
    
end